---
title: "H3"
author: "Arturo Bertero"
date: ""
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r}
#clean
rm(list = ls())

# Libraries
library(pacman)
p_load(tidyverse, here, EGAnet, psychTools, conflicted, haven, lavaan, qgraph, forcats,
       stargazer, mgm, corclass, ggcorrplot, ggpubr, nnet, ggeffects, magick, semTools,
       sjPlot, patchwork) 

#notation
options(scipen=999)

#Conflicts
conflicts_prefer(dplyr::select)
conflicts_prefer(dplyr::filter)
```

# Input

```{r}
#Load data
IPBS = readRDS((here("Input", "IPBS_imputed.rds")))

#Filter smaller dataset
att = IPBS %>% 
  select(c(abort:redis, globa, immig)) 

```

# Processing

```{r}
# Define shortnames, longnames, and groupings 
shortnames <- names(att)
longnames <- c("Abortion","Euthanasia","Same-sex marriage",
               "Redistribution", "Globalization","Immigration")

# Obj for dimensions
Totalgroup_comm <- list(
  "Cultural" = c(1:3),
  "Economic" = c(4:6)
)

Totalgroup_cols <- c("#966FD6B3","#008B8BB3")  

#remove labels
# att = sapply(att, haven::zap_labels)
```


## CCA

```{r}
#CCA 
set.seed(123)
cca_group <- cca(att, filter.significance = TRUE, filter.value = 0.05, 
                 zero.action = c("ownclass")) 

#Preliminary plots
plot_cca_1 = plot(cca_group, 1)
plot_cca_2 = plot(cca_group, 2)
plot_cca_3 = plot(cca_group, 3)
plot_cca_4 = plot(cca_group, 4)

# Add CCA membership to att and set class 2 to NA
att$cca <- cca_group$membership
att$cca[att$cca == 2] <- NA

# Same for IPBS 
IPBS$cca <- cca_group$membership
IPBS$cca[IPBS$cca == 2] <- NA
```

## SEM: do cca classes improve fit?

```{r}
#Data
SEM_data = att %>% na.omit()

# The full sample model 
model <- '
  moral =~ abort + eutha + marria
  econ  =~ redis + globa + immig
'

# Fit model to full sample (pooled)
fit_full <- sem(model, data = SEM_data, meanstructure = TRUE)

# Fit model with CCA classes (multi-group)
fit_mg <- sem(model, data = SEM_data, group = "cca", meanstructure = TRUE)

# Likelihood ratio test
anova_comp = anova(fit_full, fit_mg)

# Also check AIC/BIC directly
fitMeasures(fit_full, c("aic", "bic"))
fitMeasures(fit_mg, c("aic", "bic"))

# Save
stargazer(anova_comp, type = "html", summary = FALSE, rownames = FALSE, 
         out = here("Output", "Supplement", "Tab_5.doc"))
```


## Cor matrices

```{r}
# Data for cor matrices
att_cca_1 <- att %>% filter(cca == 1) %>% select(abort:immig)
att_cca_3 <- att %>% filter(cca == 3) %>% select(abort:immig)
att_cca_4 <- att %>% filter(cca == 4) %>% select(abort:immig)

# Correlation matrices and p-values
corr_1 <- round(cor(att_cca_1), 1)
corr_3 <- round(cor(att_cca_3), 1)
corr_4 <- round(cor(att_cca_4), 1)

p.mat_1 <- cor_pmat(att_cca_1)
p.mat_3 <- cor_pmat(att_cca_3)
p.mat_4 <- cor_pmat(att_cca_4)

# Function with 2-line title
make_corr_plot <- function(corr_matrix, title_line1, title_line2) {
  ggcorrplot(corr_matrix, type = "lower", lab = TRUE, legend.title = NULL,
             outline.col = "white", colors = c("red", "white", "blue")) +
    ggtitle(paste0(title_line1, "\n", title_line2)) +
    theme(
      plot.title = element_text(size = 14, face = "bold", hjust = 0.5, lineheight = 1.2),
      plot.margin = margin(0, 1, 0, 1, "mm"),
      axis.title = element_blank(),
      axis.text.x = element_text(size = 8, margin = margin(t = 0)),
      axis.text.y = element_text(size = 8, margin = margin(r = 0)),
      legend.position = "bottom"
    ) +
    geom_vline(xintercept = 1:ncol(corr_matrix) - 0.5, colour = "white", size = 2) +
    geom_hline(yintercept = 1:ncol(corr_matrix) - 0.5, colour = "white", size = 2)
}

# Print
cor_m_1 <- make_corr_plot(corr_1, "Class 1", "Ideologues")
cor_m_3 <- make_corr_plot(corr_3, "Class 2", "Alternatives")
cor_m_4 <- make_corr_plot(corr_4, "Class 3", "Fragmented")


# Arrange
plot_cor <- list(cor_m_1, cor_m_3, cor_m_4)
g_plot_cor <- ggarrange(plotlist = plot_cor, ncol = 3, nrow = 1, align = "hv")

# Save with a short height to force vertical expansion
ggsave(here("Output", "Supplement", "Fig_2.jpeg"),
       g_plot_cor,
       dpi = 600,
       height = 6, width = 12)

```

## Cor networks

```{r}
# Multiplot
jpeg(here("Output", "Article", "Figure_3.jpeg"), 
     height = 2480, width =3508, quality = 100)  

# Create a common layout based on the two conditions
L <- averageLayout(corr_1,corr_3,corr_4, 
                   layout = "spring")

# Matrix layout: 1 row, 3 columns
lmat <- matrix(1:3, 1, 3)

# Set layout
layout(lmat, widths = c(1, 1, 1), heights = 1)

# 1
set.seed(1)
graph1 <- qgraph(
  corr_1, layout = L, fade = TRUE,  
  labels = shortnames, vTrans = 255,
  groups = Totalgroup_comm,  
  color = Totalgroup_cols,  
  legend = FALSE, title = "1: Ideologues", title.cex = 7,
  borders = TRUE, vsize = 10, esize = 15, 
  theme = "colorblind"  
)

# 2
graph3 <- qgraph(
  corr_3, layout = L, fade = TRUE,  
  labels = shortnames, vTrans = 255,
  groups = Totalgroup_comm,  
  color = Totalgroup_cols,  
  legend = FALSE, title = "2: Alternatives", title.cex = 7,
  borders = TRUE, vsize = 10, esize = 15, 
  theme = "colorblind"  
)

# 3
graph4 <- qgraph(
  corr_4, layout = L, fade = TRUE,  
  labels = shortnames, vTrans = 255, 
  groups = Totalgroup_comm,  
  color = Totalgroup_cols,  
  legend = FALSE, title = "3: Fragmented", title.cex = 7,
  borders = TRUE, vsize = 10, esize = 15, 
  theme = "colorblind"  
)

# Finish and save the plot
dev.off()
```
```{r}
# Multiplot saved as tiff
ragg::agg_tiff(
  here("Output", "Article", "Figure_3_tiff.tiff"),
  width = 3508, height = 2480, units = "px",
  compression = "lzw", background = "white"
) 

# Create a common layout based on the two conditions
L <- averageLayout(corr_1,corr_3,corr_4, 
                   layout = "spring")

# Matrix layout: 1 row, 3 columns
lmat <- matrix(1:3, 1, 3)

# Set layout
layout(lmat, widths = c(1, 1, 1), heights = 1)

# 1
set.seed(1)
graph1 <- qgraph(
  corr_1, layout = L, fade = TRUE,  
  labels = shortnames, vTrans = 255,
  groups = Totalgroup_comm,  
  color = Totalgroup_cols,  
  legend = FALSE, title = "1: Ideologues", title.cex = 7,
  borders = TRUE, vsize = 10, esize = 15, 
  theme = "colorblind"  
)

# 2
graph3 <- qgraph(
  corr_3, layout = L, fade = TRUE,   
  labels = shortnames, vTrans = 255,  
  groups = Totalgroup_comm,  
  color = Totalgroup_cols,  
  legend = FALSE, title = "2: Alternatives", title.cex = 7,
  borders = TRUE, vsize = 10, esize = 15, 
  theme = "colorblind"  
)

# 3
graph4 <- qgraph(
  corr_4, layout = L, fade = TRUE,   
  labels = shortnames, vTrans = 255,  
  groups = Totalgroup_comm,  
  color = Totalgroup_cols,  
  legend = FALSE, title = "3: Fragmented", title.cex = 7,
  borders = TRUE, vsize = 10, esize = 15, 
  theme = "colorblind"  
)

# Finish and save the plot
dev.off()
```



## Means

```{r}
# define variable order and attitude type
var_order <- c("globa", "redis", "immig", "marria", "abort", "eutha")
att_types <- c("globa" = "Economic", "redis" = "Economic", "immig" = "Economic",
               "marria" = "Cultural", "abort" = "Cultural", "eutha" = "Cultural")

# Rescale and reshape IPBS full sample
IPBS_scaled <- IPBS %>%
  select(all_of(var_order)) %>%
  mutate(across(everything(), ~ (. - min(., na.rm = TRUE)) / (max(., na.rm = TRUE) - min(., na.rm = TRUE)))) %>%
  pivot_longer(cols = everything(), names_to = "variable", values_to = "value") %>%
  mutate(att_type = att_types[variable])

# Compute bootstrapped means and CIs for the full sample
means_full <- IPBS_scaled %>%
  group_by(variable) %>%
  summarise(
    mean = mean(value, na.rm = TRUE),
    ci_low = Hmisc::smean.cl.boot(value)[2],
    ci_high = Hmisc::smean.cl.boot(value)[3],
    att_type = first(att_type),
    .groups = "drop"
  )

# Prepare plot_data (rescaled attitudes with cca class)
plot_data <- att %>%
  mutate(cca = case_when(
    cca == 1 ~ "Class 1",
    cca == 3 ~ "Class 2",
    cca == 4 ~ "Class 3",
    TRUE ~ as.character(cca)
  )) %>%
  filter(!is.na(cca)) %>%
  mutate(cca = factor(cca, levels = c("Class 1", "Class 2", "Class 3"))) %>%
  select(all_of(var_order), cca) %>%
  mutate(across(-cca, ~ (. - min(., na.rm = TRUE)) / (max(., na.rm = TRUE) - min(., na.rm = TRUE)))) %>%
  pivot_longer(cols = -cca, names_to = "variable", values_to = "value") %>%
  mutate(
    variable = factor(variable, levels = var_order),
    att_type = att_types[variable],
    att_type = factor(att_type, levels = c("Economic", "Cultural"))
  )

# Compute bootstrapped means + Cis for each CCA class and variable
mean_ci_data <- plot_data %>%
  group_by(variable, cca, att_type) %>%
  summarise(
    mean = mean(value, na.rm = TRUE),
    ci_low = mean_cl_boot(value)$ymin,
    ci_high = mean_cl_boot(value)$ymax,
    .groups = "drop"
  )

# Colors
att_colors <- c("Economic" = "#008B8BB3", "Cultural" = "#966FD6B3")

# Relabel CCA class names
plot_data <- plot_data %>%
  mutate(
    cca = case_when(
      cca == "Class 1" ~ "Class 1\nIdeologues",
      cca == "Class 2" ~ "Class 2\nAlternatives",
      cca == "Class 3" ~ "Class 3\nFragmented",
      TRUE ~ cca
    )
  )

mean_ci_data <- mean_ci_data %>%
  mutate(
    cca = case_when(
      cca == "Class 1" ~ "Class 1\nIdeologues",
      cca == "Class 2" ~ "Class 2\nAlternatives",
      cca == "Class 3" ~ "Class 3\nFragmented",
      TRUE ~ cca
    )
  )

# Reorder variable factor 
plot_data <- plot_data %>%
  mutate(variable = factor(variable, levels = var_order))

mean_ci_data <- mean_ci_data %>%
  mutate(variable = factor(variable, levels = var_order))


#  plot
means <- ggplot(plot_data, aes(x = value, fill = att_type, color = att_type)) +
  geom_density(alpha = 0.4, size = 0.8) +
  geom_point(
    data = mean_ci_data,
    aes(x = mean, y = 0.50),  
    color = "black", size = 2, inherit.aes = FALSE
  ) +
  geom_errorbarh(
    data = mean_ci_data,
    aes(xmin = ci_low, xmax = ci_high, y = 0.50),
    height = 0.02, color = "black", inherit.aes = FALSE
  ) +
  geom_point(
    data = means_full,
    aes(x = mean, y = 0, group = NULL),
    color = "red", size = 2, inherit.aes = FALSE
  ) +
  geom_errorbarh(
    data = means_full,
    aes(xmin = ci_low, xmax = ci_high, y = 0),
    height = 0.02, color = "red", inherit.aes = FALSE
  ) +
 facet_grid(fct_rev(variable) ~ cca, switch = "y", scales = "free_y", space = "free_y")+
  scale_fill_manual(values = att_colors, name = "Attitude Type") +
  scale_color_manual(values = att_colors, name = "Attitude Type") +
  labs(
    x = "Attitude (0–1 scaled)",
    y = NULL,
    title = "Distribution of Political Attitudes by CCA Class"
  ) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    axis.title = element_text(face = "bold"),
    axis.text.x = element_text(hjust = 1),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    legend.position = "bottom"
  )

# Save
ggsave(here("Output", "Article", "Figure_4.jpeg"),
       means, dpi = 600, width = 10, height = 12)

#Save as tiff
ggsave(
  here("Output", "Article", "Figure_4_tiff.tiff"),
  means,
  dpi = 600,                
  width = 10,
  height = 12,
  device = ragg::agg_tiff,  
  compression = "lzw",      
  background = "white"      
)
```


## Boot of cca clusters

```{r}
# List of matrices to bootstrap
matrices <- list(
"att_cca_1" = att_cca_1,
"att_cca_3" = att_cca_3, 
"att_cca_4" = att_cca_4
)

# Number of bootstrap samples
n_boot <- 10000

# Function to bootstrap a matrix
bootstrap_matrix <- function(mat, n_boot) {
  lapply(1:n_boot, function(i) {
    sample_indices <- sample(1:nrow(mat), replace = TRUE)  # Sample rows with replacement
    as.matrix(mat[sample_indices, , drop = FALSE])  # Convert to matrix
  })
}

# Loop through matrices and generate bootstrap samples
boot_samples <- list()

for (name in names(matrices)) {
  cat("Bootstrapping:", name, "\n")
  boot_samples[[name]] <- bootstrap_matrix(matrices[[name]], n_boot)
}

boot_samples_cca = boot_samples
```

## Tightness of CCA 

### Cor

```{r}
# Initialize list
cor_upper_distributions <- list()

# Loop over all vote groups
for (name in names(boot_samples_cca)) {
  cat("Processing partition:", name, "\n")
  
  # Initialize vector gathering cors
  bootstrap_means <- numeric(length(boot_samples_cca[[name]]))
  
  for (i in seq_along(boot_samples_cca[[name]])) {
    data_mat <- as.matrix(boot_samples_cca[[name]][[i]])
    
    # Compute cor mat
    cor_mat <- cor(data_mat, use = "pairwise.complete.obs", method = "pearson")
    
    # Compute mean of absolute values (upper tri)
    bootstrap_means[i] <- mean(abs(cor_mat[upper.tri(cor_mat)]), na.rm = TRUE)
  }
  
  # Store all 
  cor_upper_distributions[[paste0("cor_", name)]] <- bootstrap_means
}

```

### Box plot

```{r}
# Convert to long format
cor_df <- bind_rows(
  tibble(cor = cor_upper_distributions$cor_att_cca_1, class = "Class 1\nIdeologues"),
  tibble(cor = cor_upper_distributions$cor_att_cca_3, class = "Class 2\nAlternatives"),
  tibble(cor = cor_upper_distributions$cor_att_cca_4, class = "Class 3\nFragmented")
)

#Plot
box_plot <- ggplot(cor_df, aes(x = class, y = cor)) +
  geom_boxplot(fill = "gray80", color = "black", width = 0.6, outlier.shape = NA) +
  stat_summary(fun = mean, geom = "point", size = 2, color = "black") +
  stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = 0.2, color = "black") +
  labs(
    x = "",
    y = "Mean Correlation (Bootstrap)",
    title = "Tightness by CCA Class"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5),
    axis.title = element_text(face = "bold"),
    panel.background = element_rect(fill = "white", colour = NA),
    plot.background  = element_rect(fill = "white", colour = NA)
  )

# save box_plot as figure 3_1
ggsave(here("Output", "Article", "Figure_3_1.jpeg"),
       box_plot, dpi = 600, width = 10, height = 6)

# save as tiff
ggsave(
  here("Output", "Article", "Figure_3_1_tiff.tiff"),
  box_plot,
  dpi = 600,                
  width = 10,
  height = 6,
  device = ragg::agg_tiff,  
  compression = "lzw",
  background = "white"
)

```

```{r}
#Create the final figure 3 for publication
library(magick)

# Define paths
fig1_path <- here("Output", "Article", "Figure_3_tiff.tiff")
fig2_path <- here("Output", "Article", "Figure_3_1_tiff.tiff")
output_path <- here("Output", "Article", "Figure_3_combined_tiff.tiff")

# Read images
fig1 <- image_read(fig1_path)
fig2 <- image_read(fig2_path)

# Get width of the first figure (reference)
target_width <- image_info(fig1)$width

# Resize second figure to match width (height adjusts automatically)
fig2_resized <- image_resize(fig2, paste0(target_width, "x"))

# Stack vertically
combined <- image_append(c(fig1, fig2_resized), stack = TRUE)

# Save as high-quality TIFF
image_write(
  combined,
  path = output_path,
  format = "tiff",
  compression = "lzw"
)
```


# Output

```{r}
#Descriptives on class sizes
table(IPBS$cca, useNA = "always")

#Save IPBS_with_cca
saveRDS(IPBS, here("Input", "IPBS_with_cca.rds"))

```

